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Abstract. Ground state of the dissipative two-state system is investigated by means of the Lanczos di- 
agonalization method. We adopted the Hilbert-space-reduction scheme proposed by Zhang, Jeckelmann 
and White so as to reduce the overwhelming reservoir Hilbert space to being tractable in computers. Both 
the implementation of the algorithm and the precision applied for the present system are reported in de- 
tail. We evaluate the dynamical susceptibility (resolvent) with the continued-fraction-expansion formula. 
Through analysing the resolvent over a frequency range, whose range is often called 'interesting' frequency, 
we obtain the damping rate and the oscillation frequency. Our results agree with those of a recent quantum 
Monte-Carlo study, which concludes that the critical dissipation from oscillatory to over-damped behavior 
decreases as the tunneling amplitude is strengthened. 

PACS. 75.40.Mg numerical simulation studies - 05.40.+j Fluctuation phenomena, random processes, and 
Brownian motion - 05.70.Ln Nonequilibrium thermodynamics, irreversible processes 



1 Introduction 

Effect of dissipation on quantum tunneling phenomenon 
lies out of the scope of the conventional weak coupling 
(Markovian) approximation, and has been studied exten- 
sively so far Uvm- In order to introduce dissipation, Caldeira 
and Leggett |1H] involved a thermal reservoir which con- 
sists of oscillators influencing stochastic (Brownian) fluc- 
tuations to the tunneling two-level system. Their model, 
the so-called spin-boson model, is given by the following 
Hamiltonian, 
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unit throughout this paper; u) c = 1. The transverse field 
A induces quantum coherence (tunneling amplitude) be- 
tween the spin up-and-down states, whereas the coupling 
to the reservoir disturbs the coherence. These conflicting 
effects are the central concern of the problem, which are 
apparently non-perturbative in nature. 

It is noteworthy that the spin-boson model (0) is equiv- 
alent to the anisotropic Kondo model HHSH ■ The param- 
eter a controls the strength of the longitudinal Kondo cou- 
pling, whereas A controls the transverse-coupling strength. 
Thereby, the region a < 1 (a > 1) is identified as the the 
antiferromagnetic (ferromagnetic) Kondo phase; see the 
ground-state phase diagram shown in Fig. [I]. That is, in 



i=l 



where the operators {<r Q } denote the Pauli operators and 
the operators a, and a] denote the bosonic annihilation 
and creation operators, respectively, for the i-th oscillation 
mode (i = 1 ~ N) . The set of these oscillators works as the 
above-mentioned reservoir with respect to the spin l/2er. 
The coupling coefficients {/*} should be arranged so as to 
satisfy the so-called Ohmic-dissipation condition, 
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Fig. 1. The ground-state phase diagram of the spin-boson 
model (|l|) with the Ohmic dissipation. Quantitative estimation 
of the relaxation parameters such as the damping rate and the 
oscillation frequency is the main concern of this paper. 



The dimensionless constant a controls the strength of the 
dissipation, and the cut-off frequency u> c defines the energy 



the region a < 1, the up-and-down spin states form a co- 
herent (singlet) state through a certain tunneling ampli- 
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tude, whereas in a > 1, the tunneling amplitude vanishes 
owing to the dissipation ||[l0]]. Hence, in the region a > 1, 
the ground state is degenerated doubly. The effective tun- 
neling amplitude is found to vanish at the phase boundary 
in the form fT^Q, 



A cS =A( — 



(4) 



Owing to the equivalence, we readily investigate the equi- 
librium (thermodynamic) properties by means of various 
theoretical techniques which have been developed so far 
for the Kondo problem. Note, however, the dynamical 
(non-equilibrium) properties are rather out of the scope 
of these analytical techniques; refer to the paper |H| for 
a recent analytical approach. The dynamical properties, 
especially in the frequency range ~ A e g, are the very con- 
cern of the present topic. The noninteracting-blip approxi- 
mation was invented so as to describe the relaxation 
dynamics of this particular frequency range. This approxi- 
mation is, however, justified in a rather limited parameter 
range a < 0.5. 

Hence, in order to investigate the relaxation proper- 
ties, various numerical simulations have been performed so 
far. Chakravarty and Rudnick performed quantum Monte- 
Carlo simulation p5[ ; the model they simulated is a one- 
dimensional long-range Ising model, which was derived 
||Tl|fi"7| through eliminating the reservoir (conduction elec- 
tion) degrees of freedom. They succeeded in obtaining 
the spectral function through the Pade approximation of 
the imaginary time correlation function followed by the 
analytic continuation (Wick rotation). As a result, they 
found that the long-range asymptotic form of the dynam- 
ical correlation is governed by power law. This conclusion 
was supported by Strong fig] with use of the similar nu- 
merical technique. Volker ]19[ followed the Chakravarty- 
Rudnick analysis, and reported that 'quasi-particle' pic- 
ture explains the spectral-function data. As a consequence, 
he obtained the damping rate and frequency. We utilize 
the picture to analyze our density-matrix-renormalization 
data. Costi and Kieffer |2(],|2lj used the Wilson numer- 
ical renormalization technique to calculate the spectral 
function. Their technique is particularly useful in order 
to investigate the Kondo fixed-point (very low temper- 
ature) physics definitely. Time-evolution simulation was 
performed with use of stochastic sampling (path integral) 
algorithm p^| . The method has an advantage over oth- 
ers that one can observe real-time dynamics directly. The 
stochastic-sampling error, however, grows as the evolution 
time is lengthened. 

In the present paper, for the first time, we perform 
Lanczos-diagonalization analysis of the spin-boson model 
([!]). In order to command vast assembly of the reservoir- 
oscillator modes, we adopted the density-matrix trunca- 
tion scheme proposed by Zhang, Jeckelmann and White 
p3| . (They studied the one-dimensional polaron system 
with this truncation technique.) The rest of this paper 
is organized as follows. In the next section, we propose 
an implementation of their algorithm to the spin-boson 
model (|]), and report the precision in detail. In Section ||, 



by means of this new technique, we investigate the relax- 
ation properties of the spin-boson model. We evaluate the 
dynamical susceptibility (resolvent), which is readily cal- 
culated in our scheme with use of the continucd-fraction- 
expansion formula [ p4|]25[ ] . Analyzing the analyticity of the 
resolvent, we confirm the above-mentioned 'quasi particle' 
picture [ ]19[ so as to obtain the damping rate and fre- 
quency. Our results are contrasted with the former quan- 
tum Monte-Carlo results jl9) . In the last section, we give 
summary and discussions. 

2 Application of the Hilbert-space-truncation 
algorithm to the spin-boson model 

In this section, we propose a prescription for adopting the 
Hilbert-space-reduction method to the spin-boson model. 
We then demonstrate the precision of our new scheme. 



2.1 Hilbert-space-reduction algorithm 

Even a single oscillator (boson) spans infinite-dimensional 
Hilbert space. Therefore, in order to treat an oscillator 
with the exact diagonalization method, one must truncate 
the Hilbert space inevitably. In conventional simulations, 
the boson state is represented in terms of its occupation 
number, and those states whose occupation number ex- 
ceeds a limit are disregarded (discarded). Zhang, Jeckel- 
mann and White (23j proposed an alternative representa- 
tion and a truncation criterion. Applying their scheme to 
a polaron system (one-dimensional Holstein model), they 
demonstrated that the scheme yields very precise results, 
although the dimensionality of each oscillator is reduced 
to three. This truncation is called '(numerical) renormal- 
ization. ' Their new bases are particularly efficient in those 
cases where the oscillator equilibrium position is shifted by 
a certain external force (coordinate-coordinate coupling). 
This is precisely the case of the present model @). 

In the following, we explain the detail how we adopt 
their algorithm to the spin-boson model. First, we limit 
the Hilbert-space dimensionality of each oscillator to d di- 
mensions except an oscillator with D dimensions; see Fig. 
|[ We call the d-dimensional oscillators 'small oscillators,' 
and the D-dimensional one 'big oscillator.' The choice of 
the big oscillator is to be made in sequence among the var- 
ious reservoir oscillator modes (a = 1 ~ N) . The sequence 
is continued until the relevant bases, explained below, be- 
come converged. In our experience, a few sweeps are suffi- 
cient for the convergence. The space of each small oscilla- 
tor is spanned by the above-mentioned truncated relevant 
bases, and thus the creation and the annihilation operators 
should be represented in terms of these bases correspond- 
ingly. (Before renormalization, the space is allowed to be 
of the number representation (|n),n=l~d— 1), and the 
operators are expressed in terms of this subspace.) On the 
other hand, the space of the big oscillator is spanned by 
the occupation- number representation with n = ~ D— 1. 
From this D-dimcnsional space, d relevant bases are renor- 
malized (extracted) in the following manner. 
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Fig. 2. Schematic drawing of the density-matrix- 
renormalization algorithm applied to the spin-boson model (hi). 



Finally, we mention how we chose the frequencies {uj a } 
and the coupling constants {/;}. The chose is arbitrary un- 
der the constraint eq. (0). We have uniformly distributed 
the oscillator frequency over the range < u) < lo c {= 1), 
and determined the coupling constants so as to satisfy the 
constraint. We numbered the oscillator modes (a = 1 ~ 
TV) in order of the frequency uj a (from the slowest to the 
fastest). 



2.2 Precision of the Hilbert-space-truncation algorithm 

In this subsection, we show the precision of the algo- 
rithm explained in the preceeding subsection. In Fig. |^, 
we plotted the transverse magnetization (a x ) for the sys- 
tem N — 8, A = 0.3 and a — 0.5 by means of the con- 
ventional truncation scheme; namely, the boson states are 
represented in term of the occupation number, and the 
occupation number is restricted within n < n max . (Note 



Second, we diagonalize the Hamiltonian (Q) with the 
Lanczos technique to obtain the ground state <Jo). Note 
that the Hamiltonian consists of the operators whose ma- 
trix elements arc already prepared in the above. The ground- 
state vector should be represented in the form, 



l*b> =J2^ij\i)A\j) B , 



(5) 



where the bases {|«)a} are of the big site, and the bases 
{|i)s} are °f the remaining part of the system. Thereby, 
we construct the (reduced) density matrix subjected to 
the big oscillator, 
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The new bases {uf 3 } are determined so as to diagonalize 
the density matrix, 

pvP = wpul 3 . (7) 

According to the proposal [ p6p7| , the new relevant bases 
(subspace) should be spanned by the eigenvectors up to 
the fi-th largest weight (eigenvalue) wp. That is, the new 
annihilation operator of the big oscillator is re-represented 
in terms of these d relevant eigenvectors in the following 
manner, 

N/3/3' = 'u^lA (8) 

Because the renormalization is subjected to the reduced 
density matrix, the renormalization is called 'density ma- 
trix' renormalization ||26| , |27f . 

The renormalization is continued for every oscillator 
modes (i = 1 ~ N) successively until the relevant bases 
become fixed (converged). In our experience, a few sweeps 
are sufficient for the convergence. 



Fig. 3. The transverse magnetization is plotted for the sys- 
tem with N = 8, A = 0.3 and a — 0.5 with the dimen- 
sionality T> of each oscillator varied. The plots + are evalu- 
ated with the conventional occupation-number representation. 
The occupation number is restricted within n < n m&K . Hence, 
T> — n m&K +l. The plots x are those evaluated with the density- 
matrix-truncation algorithm. The dimensionality of each small 
oscillator d is varied (X> = d) with the big oscillator dimension- 
ality fixed (D = 6). 

that the transverse magnetization indicates a degree to 
what extent the tunneling is disturbed by the reservoir 
fluctuations.) Hence, the dimensionality of each oscillator 
is given by T> — n max + 1. We observe that through in- 
creasing T>, the results converge to a certain limit. In the 
same plot, we show the results by means of the density- 
matrix- renormalization method for the same system as the 
above (JV = 8, A = 0.3 and a = 0.5). In this renormaliza- 
tion, we fixed the dimensionality of the big oscillator as 
D = 6, and varied the small-oscillator dimensionality d; 
hence, T> = d. We see that, with fewer number of bases, 
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the renormalization results converge more rapidly than 
the former occupation-number-representation results. As 
is mentioned in the previous subsection, the equilibrium 
position of each oscillator is shifted by the coordinate- 
coordinate coupling. The relevant bases of the density- 
matrix renormalization are constructed so as to take into 
account of the fluctuations from the equilibrium position, 
whereas the occupation-number bases are rather ineffi- 
cient to represent these shifted oscillator modes. 

In Fig. we show the relative error of the trans- 
verse magnetization of the density-matrix-renormalization 
method (D — 6); the error is defined as the deviation from 
the the conventional occupation-number-truncation diag- 
onalization with T> = n max + 1 = 6. The system param- 
eters are the same as those shown in Fig. y[ We observe, 
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Fig. 4. The relative error of the transverse magnetization 
with the density-matrix-renormalization method. We varied 
the number of remained bases d with D — 6 fixed. The system 
parameters are the same as those shown in Fig. |5| 

with very limited number of relevant bases (d = 2 ~ 3), 
the density-matrix renormalization reproduces the full- 
diagonalization result precisely. The relative error is sat- 
urated to ~ 10~ 7 due to the numerical round-off error 
for d > 4. The result indicates that a few relevant os- 
cillator modes are of importance. In Fig. ||, we show the 
weight Wfj of each eigenstate (i = 1 ~ D) for each reser- 
voir oscillator (a = 1 ~ N). The weight wp indicates 
the statistical weight of the state u' 3 contributing to the 
ground state. We notice that, in fact, the first few bases 
are particularly weighted, and the other states are irrel- 
evant (wp < 10~ 5 ), and to be ignored. Furthermore, it 
should be noted that these features are common to all the 
reservoir modes (i = 1 ~ N). 

To summarize, in Fig. |3|, we found that by means of 
the conventional occupation-number truncation, the result 
converges gradually as the occupation-number threshold 
n mia is increased. We see that the occupation numbers 
exceeding n m 6 are hardly excited. In the density-matrix 
renormalization, see Figs. we found that only the first 



Fig. 5. The weight vr (density-matrix eigenvalue) is plotted 
for each oscillator mode (i = 1 ~ 8). 

few states are particularly weighted. Hence, hereafter, we 
set the dimension of the big site to be D — 6, and re- 
main two relevant states (d = 2) for each oscillator mode. 
We believe that the number of the reservoir oscillators is 
prior to the number of remained bases d for each mode: 
Note that as the number of the reservoir modes is in- 
creased, the coupling constants {/,;} should be reduced 
correspondingly; the oscillators become less disturbed. In 
the following, we simulate the reservoir consisting of eigh- 
teen oscillator modes. Hence, the truncation error is ex- 
pected to be reduced furthermore from those shown in this 
subsection. 

3 Density-matrix-renormalization analysis of 
the dissipative tunneling 

With use of the method developed in the preceeding sec- 
tion, we investigate the relaxation properties in the ground 
state of the spin boson model (|j) . These properties are ex- 
tracted from the dynamical susceptibility, which is readily 
evaluated in our scheme. Our results of the relaxation co- 
efficients are contrasted with those obtained by means of 
the quantum Monte-Carlo simulation |l9[ . 

3.1 Dynamic susceptibility — 
continued-fraction-expansion formula 

In this subsection, we evaluate the dynamical susceptibil- 
ity, and compare ours with that obtained at an integrable 
point a — 0.5. The analyticity of the susceptibility is an- 
alyzed extensively in the next subsection so as to yield 
relaxation properties. 

Linear response theory states that the relaxation (non- 
equilibrium) process should be described in term of a cer- 
tain ground-state equilibrium dynamical correlation func- 
tion, unless the process is not so far from equilibrium. 
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In other words, equilibrium dynamical correlation func- 
tion does contain informations about non-equilibrium pro- 
cesses. For that purpose, we calculated the following dy- 
namical susceptibility, 



1 
2 

Im 
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(9) 
(10) 



Some might wonder that the inverse matrix of the total 
Hamiltonian appearing in eq. (^|) cannot be computed; this 
is true. The expectation value of the inverse of the Hamil- 
tonian is, however, evaluated with use of the Gagliano- 
Balseiro continued-fraction formula S.E5I, 
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Solid line: numerical result 
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Fig. 6. The spectral function S(uj) ( |l3| ) is plotted for 
A — 0.2 and a = 0.5. The solid line is our density-matrix- 
renormalization result for N = 18. The delta-function peaks 
are broadened into the Lorentz form with the width r\ = 0.022. 
The dashed line shows a rigorous result , which is available 
at the Toulouse point a = 0.5. 



where the coefficients are given by the Lanczos tri-diagonal 
elements, 



\fi+i)=H\f i )-a i \f i )-f3f\f^ 1 ), 

on = (fi\H\fi)/{fi\fi)> 

0i = (fi\fi)/(fi-x\fi-i) (#> = 0). 



(12) 



Therefore, through choosing of the Lanczos initial vector 
as l/o) = c 2 !^), we readily evaluate the dynamical sus- 
ceptibility by means of the same numerical procedure used 
in the preceeding Lanczos diagonalization. 

We expanded the formula ([Ly) up to the one-hundredth 
order. This is comparable with the iteration number needed 
to obtain the low-lying states in the Lanczos diagonaliza- 
tion. Therefore, the expansion is expected to yield precise 
result concerning low- lying frequencies. This frequency rang 
is sufficient for our purpose, because we are concerned in 
the frequency range ~ A e g. Moreover, in practice, at high 
frequencies, the magnitude of the spectral intensity is sup- 
pressed considerably. 

In Fig. ^, we plotted the spectral function, 

5H = ^M, (13) 
to 

for the parameters A = 0.2, a = 0.5 and N = 18. (The 
spectral function is related closely to the linear-response 
function.) At the point a = 0.5, the Hamiltonian ([j]) is 
reduced to being quadratic through a canonical transfor- 
mation, and the spectral intensity is calculated exactly in 
the form lEl, 



(A = (tt j '4) A 2 j '(jj c ) . The rigorous result is plotted in Fig. 
|6| as well. We observe that our numerical data reproduces 
the rigorous formula. The slight difference may be at- 
tributed to the spectral form (ps) of the reservoir around 
the cut-off frequency. We used a non-analytic reservoir 
spectrum which vanishes suddenly at the cut-off frequency. 
This difference may become irrelevant, if the tunneling 
amplitude is set to be sufficiently small compared with 
the cut-off frequency. 



3.2 Analyticity of the dynamical susceptibility 
(resolvent) and the damping properties 

In this subsection, we investigate the analyticity of the 
e dynamical susceptibility (|^), from which we extract infor- 
mations about relaxation properties. In Fig. |?j, we plotted 
the imaginary-time correlation function, 



dT(e Tn a z e- TH a z } f3 e lu ' T 



(15) 
(16) 



for the system with A = 0.2, a = 0.2 and N = 18. As 
is shown in the plot, the numerical result is fitted well by 
the 'quasi-particle' formula, 
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This fact was pointed out by Volker, who performed a 
quantum Monte-Carlo simulation; see Introduction. As is 
apparent from the definition (|9|), the (fitting) parameters 
A and u> give the damping rate and the frequency, respec- 
tively, of the 'coordinate' a z perturbed from the equilib- 
rium position. According to the formula (h_7n , the following 
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Fig. 7. The imaginary-time Green function ([l5|) is plotted 
for N — 18, A = 0.2 and a = 0.2. The solid line is our 
density-matrix-renormalization result. The plots o are those 
of the 'quasi-particle' form (^) with the damping coefficients 
A = 0.052 and ui = 0.124. 



Fig. 8. The damping coefficients A (|l^) and ui (^) are plotted; 
the system parameters are the same as those in Fig. [?]. Around 
the frequency range ~ A c fi(~ 0.13), these quantities remain 
invariant, indicating that over this range of time, the dynamics 
is described by the quasi-particle picture (fl7|). 



relations hold; 



£ (gfror 



1 m^))- 1 



(w + A) 2 



(18) 
(19) 



From these relations, we estimate the damping coefficients 
A and u>. It is one of the advantages of our scheme that 
one can differentiate the function Q(iu>), because it is ex- 
panded in the analytic form ([ll]) . in Fig. ||, we plotted the 
right hand side of eqs. ( |l~8| ) and the system param- 
eters are the same as those in Fig. 0. We see that these 
quantities are invariant actually over a certain frequency 
range ~ Z\ c ff(~ 0.13), and thus the quasi-particle picture 
( p"7| ) holds in the time range ~ A~^ . As is explained in 
Introduction, we are concerned in the physics of this par- 
ticular range of time. In consequence, we obtained the 
poles of G(oj) at lu = ±cD — iA, which are away from the 
real axis. It is noteworthy that such irreversible features 
are not transparent in the original high-order-expansion 
result ([ll]). Through the above data analysis, such the 
relaxation features become clear; we discuss this point af- 
terwards. 

In Fig. ||, we plotted the damping rate and the fre- 
quency. These damping coefficients are estimated both 
with the procedure as is shown in Fig. |^ and with the 
conventional least square fit by the function ( |l7| ) . Our re- 
sults of the density-matrix renormalization confirm the 
former quantum Monte-Carlo study |l9|j : The oscillation 
frequency uj is suppressed as the dissipation is strength- 
ened. It vanishes around a « 0.5. This point indicates 
the transition between the under-damped and the over- 
damped oscillations. By means of the bosonization tech- 
nique H , this transition point was predicted to locate at 




Fig. 9. The damping rate A and the frequency Co estimated 
by means of the density-matrix-renormalization algorithm; (a) 
A = 0.1, (b) A = 0.2 and (c) A = 0.5. These are to be con- 
trasted with those of quantum Monte-Carlo method [H. 



a = 0.5 irrespective of the tunneling amplitude strength 
A. Because in the bosonization technique, the band width 
(cut-off frequency) is supposed to be sufficiently large com- 
pared with the many-body interactions (A and a), the 
validity in the strong-coupling region is not necessarily 
assured. In Fig. ||, in fact, we see that the ui plots become 
curved convexly, as the tunneling amplitude is strength- 
ened, and accordingly it becomes evident that the transi- 
tion point locates below a = 0.5. These features were re- 
ported in the paper (l9) , and were suspected to be due to 
a systematic numerical error caused by the critical slowing 
down. In our diagonalization calculation, we are free from 
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the critical slowing down. Therefore, we conclude that for 
larger values of A, the transition point actually shifts be- 
low a — 0.5. Furthermore, we confirm the report [ fl9| that 
at a = 1/3, the damping rate and the frequency coincide. 
At this point, the peaks around u> — ±£j of spectral in- 
tensity (|l3|) merge into a single peak so that this point 
was suspected to indicate a certain phase transition. The 
present result suggests that the point is merely the situa- 
tion where the damping feature smears out the oscillatory 
behaviour, because these relative strengths exchange. 

We see that in Fig. | (a) (A = 0.1), the plots for 
a > 0.6 are rather scattered (irregular). In that region, the 
quasi-particle poles shift towards the origin of the complex 
plane (±6j — A — > 0), so that the dumping parameters 
degenerate into the slowest reservoir oscillation mode. In 
that situation, our method becomes inapplicable. Similar 
difficulty arises in the vicinity of the localization transition 
point a = 1 for larger values of A. 

We are in a position to discuss the above quasi-particle 
feature furthermore in detail. It is noteworthy that the 
form (^) implies that the analyticity of the dynamical 
susceptibility is broken along the real axis; the upper and 
lower complex planes are not continued analytically. This 
is precisely due to the reservoir effect, which drives the 
spin state to be in equilibrium, violating the time-reversal 
symmetry. In our numerical result, however, the upper 
and the lower complex planes are continued analytically, 
although along the real axis, there distribute vast num- 
ber of poles densely. This feature contradicts the above 
quasi-particle picture insisting isolated poles at ui = ±cD — 
iA. This inconsistency is simply due to the fact that our 
reservoir spectrum is not continuous, and thus the time- 
reversal symmetry is maintained. Only through the limit 
of infinite oscillator modes, these poles merge into the 
'quasi-particle' poles away from the real axis. Hence, in 
order to extract relaxation properties, one should avoid 
the real axis, along which the result is suffered signifi- 
cantly from discontinuity of the reservoir modes. Along 
the imaginary axis, as is shown in Fig. ^, the result is 
smooth and monotonic, and is much easier to be fitted 
by the quasi-particle form. That is why we examined the 
imaginary-time Green function ( p^| ) just as Volker did for 
analyzing his Monte-Carlo data computed at each Mat- 
subara frequency. We stress that our resolvent is evaluated 
(expanded) in the analytical (continued-fraction) form as 
in eq. (|TT|) . Therefore, as is shown previously in Fig. ^, one 
can obtain the spectral function (real-time Green func- 
tion) immediately through performing the Wick rotation. 
This is one of the advantages of the present scheme over 
others. 



4 Summary and discussions 

We investigated the dissipative two-state system (|l|) through 
applying the density-matrix Hilbert-space-truncation al- 
gorithm ]25| ] . An implementation of this algorithm is pro- 
posed, and the precision applied to the present model is re- 
ported in detail. We found that in our situation where the 



oscillator equilibrium position is shifted by the coordinate- 
coordinate coupling, the density-matrix renormalization 
works much better than the conventional occupation-number- 
truncation method. We have remained two relevant states 
for each oscillator mode, so that we succeeded in treating 
eighteen oscillators, keeping the truncation error consid- 
erably small. We belive that the number of tractable os- 
cillators is more crucial in observing the relaxation prop- 
erties than the 'quality' (fidelity) of each oscillator mode. 
In fact, we reproduced the dynamical susceptibility at the 
Toulouse point a — 0.5, at which a rigorous formula is 
known. We stress that, in our scheme, the susceptibility is 
evaluated in the analytical continued-fraction form ([LI]), 
which is apparently advantageous over others in perform- 
ing the analytic continuation (Wick rotation); in order to 
observe relaxation (time irreversible) properties from nu- 
merical data, which apparently possesses the time-reversal 
symmetry, we need to examine the analyticit y of the re- 
solvent (dynamical susceptibility); see Section 3.2. In fact, 
after the immediate analytic continuation to the imagi- 
nary axis, we found that the result is well fitted by the 
form (p^), confirming Volker 's finding |19 with quantum 
Monte-Carlo method. That is, the susceptibility possesses 
the so-called quasi-particle poles away from the real axis 
u> — ±w — iA. In consequence, we obtained the damping 
rate A and the frequency ui. Both agree with those of the 
quantum Monte-Carlo method |is|] . In particular, our re- 
sults confirm the former report that for larger values of 
A, the transition point between the over-damped and the 
under-damped oscillations shifts downwards (a c < 0.5). 
For studying such critical property, our method is ad- 
vantageous over the Monte-Carlo method, because our 
method is free from the critical slowing down. 
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